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We present a numerical study of the time-dependent and time-independent Gross-Pitaevskii (GP) 
equation in two space dimensions, which describes the Bose-Einstein condensate of trapped bosons 
at ultralow temperature with both attractive and repulsive interatomic interactions. Both time- 
dependent and time-independent GP equations are used to study the stationary problems. In 
addition the time-dependent approach is used to study some evolution problems of the condensate. 
Specifically, we study the evolution problem where the trap energy is suddenly changed in a stable 
preformed condensate. In this case the system oscillates with increasing amplitude and does not 
remain limited between two stable configurations. Good convergence is obtained in all cases studied. 
PACS Number(s): 02.70.Rw, 02.60.Lj, 03.75.Fi 
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I. INTRODUCTION 

Recent experiments of Bose-Einstein condensation 
(BEC) in dilute bosonic atoms (alkali and hydrogen 
atoms) employing magnetic traps at ultra-low tempera- 
tures have intensified theoretical investigations on various 
aspects of the condensate ||2|-|il|| . The properties of the 
condensate are usually described by the nonlinear mean- 
field Gross-Pitaevskii (GP) equation jl^] , which properly 
incorporates the trap potential as well as the interac- 
tion among the atoms. The GP equation in both time- 
dependent and independent forms is formally similar to 
the Schrodinger equation with a nonlinear term. The ef- 
fect of the interaction leads to the nonlinear term, which 
complicates the solution procedure. There have been sev- 
eral numerical studies of the GP equation in three space 
dimensions 

A Bose gas in lower dimensions — one and two dimen- 
sions — exhibits unusual features. For an ideal Bose gas 
BEC cannot occur in one and two space dimensions at a 
finite temperature because of thermal fluctuations ||,|l3| ■ 
The absence of BEC in one and two space dimensions 
has also been established for interacting uniform systems 
However, condensation can take place under the ac- 
tion of a trap potential |],[l4|] both for an ideal as well as 
interacting Bose gas. 

Although, there has been no experimental realization 
of BEC in two space dimensions, this is a problem of 
great theoretical and experimental interests. In a usual 
experiment of BEC in three space dimensions under the 
action of a magnetic trap the typical thermal energy 
ksT c is assumed to be much larger than energy of os- 
cillator quantum %u>, where ks is the Boltzmann con- 
stant, T c the critical temperature, and u> the oscillator 
frequency. This will allow thermal oscillation in all three 
directions. Usually, in a typical experimental situation 
the oscillator frequencies in three different directions, x, 
y, and z, are different. It is possible to obtain a quasi- 
two-dimensional BEC in a real three-dimensional trap by 



choosing the frequency in the third direction w z to satisfy 
Tilj z > ksT c > tko x , huiy. In that case the energy for ther- 
mal fluctuation is much smaller than the oscillator energy 
in the z direction. Consequently, any motion in the z di- 
rection will be frozen and this will lead to a realization 
of BEC in two space dimensions. The main features of 
BEC in two dimensions under the action of a harmonic 
trap has been discussed by Mullin recently Hj. Also, 
there has been consideration of BEC in low-dimensional 
systems for particles confined by gravitational field or by 
a rotational container p5| . Possible experimental con- 
figurations for BEC in spin-polarized hydrogen in two 
dimensions are currently being discussed |||| . 

Recent numerical studies of the GP equation in three 
space dimensions [Q-Q] in time-independent and time- 
dependent forms have emphasized that extensive care 
in numerical integration is needed to obtain good con- 
vergence. With the viability of experimental detec- 
tion of BEC in two space dimensions ||, here we per- 
form a numerical study of the time-dependent and time- 
independent GP equation in two space dimensions for an 
interacting Bose gas under the action of a harmonic oscil- 
lator trap potential. The interatomic interaction is taken 
to be both attractive and repulsive in nature. 

The nonlinear time-dependent and time-independent 
GP equations can be compared with the corresponding 
two types of the linear Schrodinger equation. The sta- 
tionary states in both cases have a trivial time depen- 
dence of the form ^(r, t) = exp(—iEt/K)^(r) where E is 
the parametric energy and t the time. As is well known 
the time-independent form of these equations determines 
the stationary function 'l'(r), as in the hydrogen-atom 
problem. The time-dependent Schrodinger equation can 
also be directly solved to obtain the full time-dependent 
solution in the case of the stationary problems, from 
which the trivial time dependence exp(—iEt/h) can be 
separated. In fact, the time-dependent methods have 
been successfully used for the bound-state calculation in 
many areas of computational quantum chemistry p6[ . 
This way of extracting the stationary solution from the 
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linear time-dependent Schrodinger equation continues 
as a powerful technique in the case of nonlinear time- 
dependent GP equations. 

In this paper we solve the stationary BEC problem 
in two dimensions using both the time-dependent and 
time-independent GP equations in the cases of attrac- 
tive and repulsive interatomic interactions and compare 
the two types of solutions. The time-independent GP 
equation is solved by integrating it with the Runge- 
Kutta rule complimented by the known boundary con- 
ditions at origin and infinity p7[ . The time-dependent 
GP equation is solved by discretization and Gauss elimi- 
nation method with the Crank- Nicholson-type rule com- 
plimented again by the known boundary conditions 
|p7| . We find that both the time-dependent and time- 
independent approaches lead to good convergence for the 
stationary bound-state problem of the condensate. We 
also compare these solutions with the Thomas-Fermi ap- 
proximation in the case of repulsive interatomic interac- 
tion. 

In addition to obtaining the solution of the stationary 
problem the time-dependent GP equation can be used to 
study the intrinsic time-evolution problems with nontriv- 
ial time dependence and in this paper the time-dependent 
approach is also used to study some evolution problems. 
Specifically, we study the effect of suddenly altering the 
trapping energy on a preformed condensate. We find 
that in this case instead of executing sinusoidal oscilla- 
tions between the stable initial and final configurations as 
in standard time-evolution problems governed by the lin- 
ear Schrodinger equation, the condensate executes oscil- 
lations around the stable initial and final configurations 
with ever-growing amplitude. 

In Sec. II we describe the time-dependent and time- 
independent forms of the GP equation. In Sec. Ill we 
describe the numerical method in some detail. In Sec. 
IV we report the numerical results and finally, in Sec. V 
we give a summary of our investigation. 



II. NONLINEAR GROSS-PITAEVSKII 
EQUATION 

At zero temperature, the time-dependent Bose- 
Einstein condensate wave function W(r, r) at position r 
and time r may be described by the self-consistent mean- 
field nonlinear GP equation In the presence of a 
magnetic trap this equation is written as 



OT 



*(r,r) = 0. 



(2.1) 



Here m is the mass of a single bosonic atom, N the num- 
ber of atoms in the condensate, mw 2 r 2 /2 the attractive 
harmonic-oscillator trap potential, lj the oscillator fre- 
quency, and g the strength of interatomic interaction. A 



positive g corresponds to a repulsive interaction and a 
negative g to an attractive interaction. The normaliza- 
tion condition of the wave function is 



dr|*(r,i)| s 



1. 



(2.2) 



For a stationary solution the time dependence of the 
wave function is given by ^(r, r) = exp(— i/i7"/7i)\&(r) 
where \x is the chemical potential of the co nde nsate. If we 
use this form of the wave function in Eq. (2J), we obtain 
the following stationary nonlinear time-independent GP 
equation |12|: 



h 1 

V 2 + -mu 2 r 2 + gNl^ir^* 2 

2m 2 



)| 2 -M 



tf(r) = 0. (2.3) 



The time-dependent equation ( |2.l| ) is equally useful for 
obtaining a stationary solution with trivial time depen- 
dence as well as for studying evolution processes with 
explicit time dependence. 

Here we shall be interested in the spherically symmct- 
ric s olutio n ^( r, r) = ip(r,r) — (p(r) exp(— i\ir jti) toEqs. 
(|2.lD and (2.3), which can be written, respectively, as 
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gN\^(r,t)\ 2 
ip(r,r) = 0, 



(2.4) 



Tx 1 d d 1 2 2 i\ t i / \i2 

2^rTr r d-r + 2 mUJr 



ip(r) = 0. 



(2.5) 



The above limitation to the spherically symmetric so- 
lution (in zero angular momentum state) reduces the 
GP equations in two physical space dimensions to one- 
dimensional differential equations. We shall study nu- 
merically these one-dimensional equations in this paper. 

As in Ref . Q , it is convenient to use dimensionless vari- 
ables defined by x = r/a n0 , and t = tw/2. where «j KJ = 

\JUj (raw), a = /j,/(%uj), ip(x) — a\ l0 \/2mgNif{r) /%, and 
i Mx, t) = au Q y / 27ry(r, r). In terms of these variables Eqs. 
(|2.4|) and (2.5) becomes, respectively, 
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--— x— +x +cn\ip(x,t)\ — 
x ox ox at 



rP(x,t) = 0, (2.6) 



—x— + x 2 + c\ip(x)r - 2a 

x ax ax 



ijj(x) = 0. (2.7) 



where n = mgN j (irk 2 ) is the reduced number of particles 
and c = ±1 carries the sign of g: c — 1 corresponds to 
a repulsive interaction and c = — 1 corresponds to an 
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attractive interaction. The normalization condition (2.2 
of the wave functions become 



\tp(x,t)\ 2 xdx 



1 



\ip{x)\ 2 xdx. (2. 



We shall be using these two slightly different normaliza- 
tions of the time-dependent and time-independent wave 
functions for future numerical convenience. 

An interesting property of the condensate wave func- 
tion is its mean-square radius defined by 



1 f°° 

x 2 \ip(x, t)\ 2 xdx = — I x\ip(x)\xdx. 
i nj 



(2.9) 



III. NUMERICAL METHOD 

A. Boundary Condition 

Both in time-dependent and time-independent ap- 
proaches we need the boundary conditions of the wave 
function as x — ► and oo. For a confined condensate, 
for a sufficiently large x, ip(x) must vanish asymptoti- 
cally. Hence the nonlinear term proportional to |V>(x)| 3 
can eventually be neglected in the GP equation for large 
x and Eq. (2.7) becomes 



1 d d 2 

x dx dx 



2a 



ijj{x) = 0. 



(3.1) 



This is the equation for the oscillator in two space di- 
mensions in the spherically symmetric state with solu- 
tions for a = 1, 3, 5, ... etc. A general classification of all 
the states of such an oscillator is well under control |L8) . 
In the present BEC problem, Eq. (3.1) determines only 
the asymptotic behavior. If we consider Eq. ( |3.1| ) as a 
mathematical equation valid for all a and large x, the 
asymptotic form of the physically acceptable solution is 
given by 



lim ip{x) = Nq exp 



-— + (a — 1) In a; 



(3.2) 



where Nc is a normalization constant. Equation (3.2) 
leads to the following asymptotic log-derivative 



lim 



a — 1 



(3.3) 



which is independent of the constant Nq and where the 
prime denotes derivative with respect to x. 

Next we consider Eq. (2.7) as x — > 0. The nonlinear 



term approaches a constant in this limit because of the 
regularity of the wave function at x = 0. Then one has 
the following usual conditions 



ip(0) = constant, 



V'(0) - 0, 



(3.4) 



as in the case of the harmonic oscillator problem in two 
space dimensions Both the small- and large- x be- 

haviors of the wave function will be necessary for a nu- 
merical solution of the GP equation in time-dependent 
and time-independent forms. 



B. Time-Dependent Approach: Evolution and 
Stationary Problems 

First we describe the nu meri cal method for solving the 
time-dependent equation (2.6). For a numerical solu- 
tion it is convenient to make the substitution ip(x,t) = 
cj)(x,t)/x in this equation, when this equation becomes 



dx 2 
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x dx 
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4>{x,t) = o. 



(3.5) 



A convenient way to solve Eq. (3J3) numerically is to 
discretize it in both space and time and reduce it to a 
set of algebraic equations which could then be solved by 
using the known asymptotic boundary conditions. We 
discretize this equation by using a space step h and time 
step A with a finite difference scheme using the unknown 
which will be approximation of the exact solution 
<fi(xj,tk) where Xj — jh and tk = kA. As Eq. (3.5) in- 
volves both time and space variables it can be discretized 
in more than one way. The time derivative in Eq. ( |3.5| ) 
involves the wave function at times tk and tk + A. As A 
is small, the time-independent operations in this equa- 
tion can be discretized by using the wave-function com- 
ponents at time tk or tfc+i = tk + A. If one uses the 
wave- function components at time tk, Eq. ( |3.5[ ) is dis- 
cretized as 
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(3.6) 



This is an explicit differencing scheme, since, given <f> at 
tk it is straightforward to solve for cf> at tk+i ]17| |. One 
should start with an approximately known solution at 
tk and propagate it in time until a converged solution 
is reached. We confirm in our study that this simple 
scheme leads to slow convergence and large unphysical 
oscillations in the solution. 

One c an e xpress the derivatives on the right-hand-side 
of Eq. (3.6) in terms of the variables at time tk+i p7| . 
Then the unknown 4>^ +1 appears on both sides of the 
equation and one has an implicit scheme. We find that 
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the implicit scheme improves substantially the numerical 
accuracy and convergence rate. However, we find after 
s om e experimentation that if the right-hand-side of Eq. 
(3.6) is averaged over times tk and tk+i one has the best 
convergence. This is a semi-implicit scheme based on 
the Crank-Nicholson scheme for discretization jl9). We 
use the fol low ing rule to discretize the partial differential 
equation (|J) 0® 



A 



1 

'2ft 5 
1 



2^ fc+1 + ^+ 1 ) 



1 



mi 



cn- 



(3.7) 

Similar discretization rule has been used for the solution 
of the GP equation in three space dimensions j|. The 
first and second space derivatives of the wave function as 
well as the the wave function itself have been approxi- 
mated by the average over their values at the initial time 
tk and the final time ffc+i. This procedure leads to ac- 
curate and stable numerical results. Consid erin g that 
the wave function is known at time tk, Eq. (3.7) is an 
equation in three unknowns — ^j+i > and ytt\- In 
a lattice of N points Eq. (3.7) represents a tridiagonal 
set for j — 2,3, (N — 1). This set has a unique so- 
lution if the wave functions at the two end points <p\ +1 



and 



,fc+i 
"n 



are known. In the present problem these values 



at the end points are provided by the known asymptotic 
conditions. The tridiagonal set of equations is solved 
by the Gauss elimination method and back substitution 
p7[ using a typical space step h = 0.0001 and time step 
A = 0.03. Although, the iterative method should work 
for any value of A, we found the convergence to be faster 
with this value of A and we used this value throughout 
the present investigation. 

The time-dependent method could be used to study 
stationary as well as time-evolution problems. First we 
consider the stationary problem. For the ground and 
the first excited states of the condensate we start with 
the following analytically kn own wave functions of the 
harmonic oscillator problem (|3.l|) |18| : 



4>{x) — xip(x) = v / 2a;exp(— x 2 /2), 
<j){x) = xnp(x) = V2x(l - x 2 ) exp(-x 2 /2), 



(3.8) 



(3.9) 



respectively, at an initial time t = 0. We then repeat- 
edly propagate these solutions in time using the Crank- 
Nicholson- type algorithm (3.7). The boundary condition 
(|3.4[), that 0(0) = 0, is implemented at each time step 



fll7| . Also, the solution at each time step will satisfy the 
asymptotic condition ( |3.2|) . Starting with cn = 0, at 
each time step we increase or decrease the nonlinear con- 
stant cn by an amount Ai typically around 0.01. This 
procedure is continued until the desired final value of cn 
is reached. Then the final solution is iterated several 
times (between 10 to 40 times) to obtain a stable con- 
verged result. The resulting solution is the ground state 
of the condensate corresponding to the specific nonlinear 
constant cn. We found the convergence to be fast for 
small \cn\. However, the final convergence of the scheme 
breaks down if \cn\ is too large. In practice these difficul- 
ties start for cn > 20 for the ground state for a positive 
c (repulsive interaction) in a computational analysis in 
double precision. For an attractive interaction there is 
no such problem as the GP equation does not sustain 
a large nonlinearity \cn\ as we comment in detail in the 
next section. 

As the time dependence of these stationary states is 
trivial — t(i(x,t) — ip(x) exp(— i2ott) — the chemical po- 
tential a can be obtained from the propagation of the 
converged ground-state solution at two successive times, 
e.g., ip(x,tk) and ip(x, ffc+i). From the numerically ob- 
tained ratio ip(x,tk)/ip(x,tk+i) = exp(i2oA) a can be 
obtained as the time step A is known. 

The time-dependent method could also be used to 
study evolution problems. One such evolution problem 
describes the fate of the condensate if the trap poten- 
tial is removed or altered suddenly after the formation of 
the condensate. As a stable condensate is formed under 
the action of the trap potential, after a sudden change in 
the trap potential, the condensate will gradually modify 
with time. To study the time evolution of a condensate 
wave function as the trap is removed or altered suddenly, 
we have to start the time evolution of the known precal- 
culated wave function of the condensate with the initial 
trap potential and allow it to evolve in time using the 
time-dependent GP equation with the full nonlinearity 
but with the altered trap potential, which could be zero. 



C. Time-Independent Approach 



The time-independent GP equation (2.7) has the fol- 
lowing structure 



y = G(x,ip(y)), 



(3.10) 



with y — xip' ', where the prime denotes the x deriva- 
tive. W ith this realization, a numerical integration of 
Eq. (^[t]) can be implemented using the following four- 
point Runge-Kutta rule [ pTpC| ] in steps of h from Xj to 

Xj+l 



ipj+2 = ipj+i + hip' j+1 , 
1 



Zj+iipj+i = Xjipj + -(s + 2si + 2s 2 + S3), 
where 



(3.11) 
(3.12) 
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s = hG(xj,ipj), 



si = hG 
s 2 = hG 
s 3 = hG 



h , h(xM + a /2) 

1 



2(xj + h/2) 
Hxjtp'j + si/2) 



x i + 2>^ + 2{ Xj + h/2) 



h{xM + s 2 ) 



(xj + h) 



(3.13) 
(3.14) 

(3.15) 

(3.16) 



Equation (2.7) is integrated numerically for a given a us- 
ing this algorithm starting a t th e origin (x — 0) with the 
initial boundary condition ( |3.4| ) with a trial ^>(0) and a 
typical space step h = 0.0001. The integration is propa- 
gated to x = imax, where the asymptotic condition (3.2) 
is valid. The agreement between the numerically calcu- 
lated log-de riva tive of the wave function and the theoret- 
ical result (^3) was enforced to five significant figures. 
The maximum value of x, up to which we needed to in- 



tegrate (2.7) numerically for obtaining this precision, is 
Xm&x = 5. If for a trial i/j(0), the agreement of the log- 
derivative can not be obtained, a new value of ^(0) is to 
be chosen. The proper choice of t/)(0) was implemented 
by the secant method. Even with this method, some- 
times it is difficult to obtain the proper value of ip(Q) for 
a given a. Unless the initial guess is "right" and one is 
sufficiently near the desired solution, the method could 
fail, specially, for large \cn\ and lead numerically to either 
the trivial solution ip{ x ) — or an exponentially diver- 
gent nonnormalizable solution in the asymptotic region. 



IV. NUMERICAL RESULT 



A. Stationary Problem 



First we consider the ground-state solution of Eq. (2.7) 
for different a in cases of both attractive and repulsive 
interactions using the time-independent method. In the 
presence of the nonlinearity, for attractive (repulsive) in- 
teratomic interaction, the solutions of the GP equation 
for the ground state appear for values of chemical poten- 
tial a < 1 (a > 1). The relevant parameters for the so- 
lutions — the wave-function at the origin ^(0), reduced 
number n, and mean-square radii (a; 2 ) — are listed in 
Table I. The numerical integration was performed up to 
a^max = 5 wit h h = 0.0001 where the asymptotic bound- 
ary condition (3J5) is implemented. 

Using the known tabulated values of n in each case 
we also solved the time-dependent GP equation and 
the wave functions and energies so calculated agree well 
with the respective quantities calculated with the time- 
independent approach. The solutions were obtained us- 
ing space step h — 0.0001, time step A = 0.03 and the 
parameter Ai ft! 0.01. The largest value of x used in 
discretization (3.7) is Xmax = 10. The wave functions 
for different values of a (and n) for the attractive and 
repulsive interparticle interactions for the cases shown in 



Table I are exhibited in Figs. 1(a) and 1(b), respectively, 
where we plot ip( x ) versus x using the time-dependent 
and time-independent approaches. The curves in Figs. 
1(a) and 1(b) appear in the same order as the rows in 
Table I and it is easy to identify the corresponding val- 
ues of a from the values of ip(0) °f each curve. From 
Figs. 1(a) and (b) we find that the nature of the wave 
function for these two cases are quite different. However, 
the wave functions calcula^d^ with time-dependent and 
time-independent approaches agree reasonably with each 
other. 
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Fig. 1. Ground-state condensate wave function ip(x) 
versus x for (a) attractive and (b) repulsive interparticle 
interactions using the time-dependent (dashed line) and 
time-independent (full line) approaches. The parameters 
for these cases are given in Table I. In the time-dependent 
method we used time step A = 0.03, Ai = 0.01, space 
step h = 0.0001 and Xmax = 8, in the time-independent 
method we used space step h — 0.0001 and Xmax = 5. In 
the case of the repulsive interparticle interaction we also 



show the solution (4.1) corresponding to the Thomas- 
Fermi approximation (dashed-dotted line). The curves 
appear in same order as in Table I. with the lowermost 
curve corresponding to the first row. 

In the absence of previous solutions of this problem 
we compare the stationary solutions in the repulsive case 
(c = 1) with those obtained via a well-known approxi- 
mation, e.g., the Thomas- Fermi approximation. In t his 
approximation the kinetic energy term in Eq. (2.7) is 



neglected and one has the following simple approximate 
solution 
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ip(x) — \J 2a - 



(4.1) 



for x 2 < 2a and zero otherwise. In Fig. 1 (b) we also plot 



the Thomas- Fermi approximation (4.1). We find that as 
expected, for a large condensate, this approximation is a 
reasonable approximation. However, it turns out to be a 
bad approximation for a small condensate. 

Table I: Pa ram eters for the numerical solution of the 
GP equation ( ^2 . T[) for c = ±1 for the ground state wave 
function. The first four columns refer to the attractive 
interaction c — — 1 and the last four columns refer to the 
repulsive interaction c = 1. 
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0.8 


0.9185 


0.3663 


0.9030 


1.2 


0.8719 


0.4353 


1.1027 


0.4 


1.6795 


0.9147 


0.7297 


1.6 


1.4415 


1.5276 


1.3219 


-0.4 


2.8255 


1.4798 


0.4757 


2.2 


1.9276 


3.7509 


1.6741 


-2.0 


4.6249 


1.7695 


0.2400 


3.0 


2.3626 


7.8377 


2.1679 


-4.0 


6.3252 


1.8319 


0.1385 


4.0 


2.7786 


14.7609 


2.8041 



It is appropriate to comment on the numerical ac- 
curacy of the present time-dependent and independent 
methods, which seems to be limited typically by the dif- 
ference between the time-dependent and independent so- 
lutions in Fig. 1. When the solution can be obtained 
numerically, as in the cases shown in Table I, the time- 
independent method can yield very accurate results. This 
accuracy can be increased by controlling the space step 
h and x m ax- This is not so in the case of the time- 
dependent method, where the numerical result exhibits 
small periodic oscillation after iteration specially for large 
values of \cn\ which we detail below. 

The numerical solution of the time-dependent method 
is independent of the space step h provided that a typical 
value around h = 0.0001 is employed as in the present 
study. No visible difference in the solution is found if h is 
increased by a factor of 2 or 3. However, the solution is 
more sensitive to the number of time iterations, specially, 
for a large value of \cn\, for a fixed integration time step 
A or the step Ai by which the nonlinear constant in the 
GP equation is increased at each time step until the final 
value of cn is reached. We show this variation in Figs. 
2 (a) and (b) where we plot \ip(x, T)\ as a function of 
reduced time T = i/0.03 for x = and 2 for different 
choices of A and Ai in the repulsive case for n = 3.7509 
and a = 2.2 corresponding to the fourth row of Table I. 
The zero of reduced time T is made to coincide with the 
time step tk at which the full nonlinear constant cn is 
obtained for the first time during iteration. This choice 
of time will allow us to compare the fluctuations of the 
solution during the time propagation of the full GP equa- 
tion. In Fig. 2(a) we present our results for A = 0.03 
and for Ai = 0.018754 and 0.0046886. In Fig. 2(b) we 
present our results for Ai = 0.0046886 and for A = 0.03, 
and 0.05. From Figs. 2 (a) and (b) we find that there is 
numerical oscillation of the solution with time in this ap- 
proach which is independent of small variations of A near 



0.03 and Ai around 0.01. These oscillations determine 
the numerical error of the time-dependent approach and 
become larger when we employ a A very different from 
0.03, or Ai very different from 0.01. The oscillations 
can really be large if an improper value of step A or Ai 
is choosen as can be seen from Fig, 2(b) for A = 0.05. 
The results remain stable if we reduce these steps up to 
A ~ 0.01 and Ai ~ 0.003. For very small Ai and A 
accumulative errors also increase. This accumulative nu- 
merical error increases as the number of iterations is very 
large (several thousands) and a large number of iterations 
is needed to cover a given time interval with a small time 
step A. 

Table II: Amplitude of oscillation A(x, T) (in units of 
0.01) of \i/j(x,T)\ at different times T for x = and 2 
calculated with A = 0.03 and A x = 0.0046886 in the 
repulsive case for cn = 3.7509. The average value of 
converged |^(0,T)| = 1.9310 and |V>(2,T)| = 0.6667. 
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\A(0,T)\ 
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1.27 


294 
3.32 
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7.35 
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\A(2,T)\ 




0.77 


162 
0.95 


291 
1.13 


400 
1.33 
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Fig. 2. 

\ip(x,T)\ versus reduced time T = i/0.03 for x = 
and 2 in the repulsive case for the nonlinear constant 
cn = 3.7509 for (a) A = 0.03, and A t = 0.0046886 (full 
line), and 0.018754 (dashed line) and (b) Ax = 0.0046886 
and for A = 0.03 (full line), and 0.05 (dashed line). The 
zero of T is taken to be the time at which the full non- 
linearity is achieved for the first time. 
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We show a quantitative account of the above oscillation 
in Table II where we plot the maximum error in \ip(x, T)\ 
(amplitude of oscillation of \ip(x,T)\) for x — and 2 at 
different times calculated with steps A = 0.03 and Ai = 
0.0046886. We find that the error increases slowly, but 
not necessarily monotonically, with time. The average 
value of the converged |-0(O,T)| is 1.9310 and that for 
|^(2, T)\ is 0.6667. The maximum deviations from these 
values as shown in Table II do not occur at the same 
values of T. We find from Table II that for small T(~ 0) 
the maximum average error in \ip(x,T)\ is about 1%. For 
T ~ 800 this maximum average error could be as high as 
4%. As these errors are oscillating with time, at a given 
T this error could be smaller or even zero. Considering 
that we are dealing with nonlinear equations these errors 
are well within the acceptable limits. The errors shown in 
Table II would also be the typical errors in time-evolution 
problems which we study in the next subsection. 

For repulsive interaction, it was increasingly difficult 
to find the solution of the GP equation using both time- 
dependent and time-independent methods for larger non- 
linearity than those reported in Figs. 1 (a) and (b) . The 
inputs of the time-independent method are a and an ap- 
propriate ip(0). In this method it became difficult (or 
impossible) to find the appropriate ijj(0) and find a solu- 
tion for large cn(> 20). For large nonlinearity the secant 
method led to radially excited state for the appropriate 
"0(0). In the time-dependent method the only input is the 
value of cn. For a large cn in the repulsive case, the nu- 
merically obtained solution for the wave function shows 
many oscillations and is clearly unacceptable physically. 
A Crank-Nicholson-type approach was also used to solve 
the GP equation in three space dimensions The nu- 
merical instability also set a limit in that investigation in 
finding stationary ground-state solution for large values 
of nonlinearity. 

For attractive interparticle interaction, the wave func- 
tion is more sharply peaked at x = than in the case 
of the repulsive interparticle interaction and one has a 
smaller reduced number n and mean square radius (x 2 ). 
In this case we find from Table I that with a reduction of 
the chemical potential a, the reduced number n increases 
slowly and the mean square radius (x 2 ) decreases rapidly, 
so that the density of the condensate p = n/(x 2 } tends 
to diverge as n tends to a maximum value nmax- The 
increase in density lowers the interaction energy. The 
kinetic energy of the system is responsible for the sta- 
bilization. As the central density increases further for 
stronger attractive interparticle interaction, kinetic en- 
ergy can no longer maintain equilibrium of the system 
and the system collapses. Consequently, for n > n m ax, 
there is no stable solution of the GP equation. Numeri- 
cally, from a plot of n versus 1/p we find this maximum 
number consistent with p _1 = to be 



nmax = r) N n 



1. 



(4.2) 



the condensate increases in size as the number of par- 
ticles in the condensate increases. These behaviors of 
the Bose-Einstein condensate in two dimensions were 
also noted in three dimensions jnj[n]. However, in 
three dimensions the correspon ding maximum value was 
nmax = 4iV max |a|/a ho « 2.30 g(j. 

Both the time-dependent and time-independent ap- 
proaches are equally applicable for spherically-symmetric 
radially excited states. For the first excited state, with 
one node in the wave function, we verified that the con- 
vergence was as good as in the ground-state case reported 
here. However, it is a routine study and we do not report 
the results here. 
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Fig. 3. Condensate wave function ip(x) in the repul- 
sive case at different times T — t/0.03 for an expanding 
condensate after the trap is removed suddenly at T = 0. 
The initial condensate has n = 8, and the time evolution 
is performed using time step A = 0.03, Ai = 0.01, space 
step h = 0.0001 and £max' B = JL5. 




There is no such limit on n in the repulsive case. In 
that case with the increase of the chemical potential a 



Fig. 4. The central probability density |?/>(0)| 2 in 
the repulsive case at different times T = i/0.03 for an 
oscillating condensate with n = 8 after the trap energy 
is suddenly reduced to half at T = 0. The \ip{Q)\ 2 for 
condensates corresponding to the initial and final traps 
are denoted by the two straight lines. The time evolution 
is performed using time step A = 0.03, Ai = 0.01, space 
step h = 0.0001, and x m &x = 15. 
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Next we consider two time-evolution problems using 
the time-dependent method. We consider the ground 
state in the repulsive case with cn = 8. In the first prob- 
lem, at t = 0, the trap is suddenly removed. In the 
second, at t — 0, the trap energy is suddenly reduced to 
half of the starting value. In both cases we study how the 
system evolves with time by solving the time-dependent 
GP equation using time step A = 0.03, Ai ~ 0.01 and 
space step h = 0.0001. Both these problems are intrinsic 
time-dependent problems and can be studied numerically 
and experimentally. 

The condensate cannot exist in the absence of the trap. 
In the first case after the trap is removed at t = 0, the 
radius of the condensate increases and the wave function 
extends over a larger region of space. We solve the time- 
dependent GP equation at different times. In Fig. 3 we 
plot the wave function at different reduced times T = 
i/0.03. The condensate increases in size monotonically 
with time and eventually F Sisappears. 
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Fig. 5. Condensate wave function i/j(x) at different 
times T = t/0.03 for an oscillating condensate after the 
trap energy is suddenly reduced to half at T — 0. All 
parameters are the same as in Fig. 3. 

In the second problem at t = 0, we reduce the trap 
energy suddenly to half of the initial value correspond- 
ing to a stable final configuration for the condensate in 
the repulsive case. The system is now found to oscillate 
between the initial and final stationary states. In the ab- 
sence of the nonlinearity, the system executes sinusoidal 
oscillations between the two stable configurations. How- 
ever, in this nonlinear problem the system executes os- 
cillations with evergrowing amplitude. To illustrate this 
oscillation we plot in Fig. 4 the the central probabil- 
ity density \ip(0)\ 2 versus reduced time T = t/0.03. The 
|^(0)| 2 for condensates corresponding to the initial and 
final trap energies are denoted by the two straight lines. 
We see that the oscillation increases with time. In our 
numerical study we find that after a very large number of 
iterations (several thousands) the amplitude may become 
very large. However, we do not know if this result makes 
sense physically as the cumulative numerical error of the 
type shown in Fig. 2 will also grow after a very large 
number of iterations, which will possibly invalidate our 



conclusion. However, the solution presented in Fig. 4 is 
stable numerically and is the acceptable physical solution 
of the problem after a small number of iterations. This 
interesting behavior can possibly be observed experimen- 
tally and deserves further theoretical and numerical stud- 
ies. In Fig. 5 we plot the wave functions of the system at 
different times which have very acceptable and smooth 
behavior. As the number of particles of the system con- 
tinues fixed, the wave functions of smaller amplit ude s 
have larger spacial extension [mean squa re radius (|2.S| )] 
so that the normalization condition (2.8) is preserved. 



V. SUMMARY 

In this paper we present a numerical study of the 
Gross-Pitaevskii equation for BEC in two space dimen- 
sions under the action of a harmonic oscillator trap po- 
tential for bosonic atoms with attractive and repulsive 
interparticle interactions using time-dependent and time- 
independent approaches. Both approaches are used for 
the study of the stationary problem. In addition some 
evolution problems were studied by the time-d epen dent 
appr oach. We derive the boundary conditions ( |3.3|) and 
(3.4) of th e so lution of the dimensionless GP equations 
(2.6) and (2/7). These boundary conditions are used for 
the solution of the stationary problem using both the 
time-dependent and time-independent approaches. 

The time-dependent GP equation is solved by dis- 
cretizing it using a Crank-Nicholson-type scheme, 
whereas the time-independent GP equation is solved by 
numerical integration using the four-point Runge-Kutta 
rule. In both cases numerical difficulty appears for large 
nonlinearity (cn > 20). For medium nonlinearity, the ac- 
curacy of the time-independent method can be increased 
by reducing the space step h. However, the solution of 
the time-dependent approach exhibits intrinsic oscilla- 
tion with time iteration which is independent of space 
and time steps used in discretization. 

The ground-state wave function is found to be sharply 
peaked near the origin for attractive interatomic inter- 
action. For a repulsive interatomic interaction the wave 
function extends over a larger region of space. In the 
case of an attractive potential, the mean square radius 
decreases with an increase of the number of particles in 
the condensate. Consequently, a stable solution of the 
GP equation can be obtained for a maximum number of 
particles in the condensate as given in Eq. 



(4.2) 



In addition to the stationary problem we studied two 
evolution problems using the time-dependent approach. 
A stable bound state is considered and the trap poten- 
tial is suddenly removed or reduced to half at t = 0. 
If the trap is removed suddenly, the system gradually 
and monotonically increases in size with time and even- 
tually it disappears occupying the whole space with zero 
density. If the trap energy is suddenly reduced to half, 
the system oscillates around the two stationary positions. 
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The amplitude of the oscillation continue to increase with 
time. This behavior is interesting and can be studied ex- 
perimentally in the future. 

The work is supported in part by the Consclho 
Nacional de Desenvolvimento Cicntffico e Tecnoland 
Fundagao de Amparo a Pcsquisa do Estado de Sao Paulo 
of Brazil. 
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